home *** CD-ROM | disk | FTP | other *** search
/ Pascal Super Library / Pascal Super Library (CW International)(1997).bin / MATH / MATH1 / LEAST1.PAS < prev    next >
Pascal/Delphi Source File  |  1985-04-03  |  4KB  |  182 lines
  1. {$S+ }
  2. program least1;        { --> 191 }
  3. { Pascal Program to perform a liner least-squares fit using a parabolic }
  4. {    curve. Aeperate procedure PLOT needed }
  5.  
  6. const    maxr    = 20;
  7.     maxc    = 3;
  8.  
  9. type    ary    = array[1..maxr] of real;
  10.     arys    = array[1..maxc] of real;
  11.     ary2s    = array[1..maxc,1..maxc] of real;
  12.  
  13. var    x,y,y_calc    : ary;
  14.     coef        : arys;
  15.     nrow,ncol    : integer;
  16.     correl_coef    : real;
  17.  
  18.  
  19. external procedure cls;
  20.  
  21. procedure get_data(var x,y: ary;
  22.            var nrow: integer);
  23. { get values for n and arrays x,y }
  24.  
  25. var    i    : integer;
  26.  
  27. begin
  28.   nrow:=9;
  29.   writeln;
  30.   for i:=1 to nrow do x[i]:=i;
  31.   y[1]:=2.07;    y[2]:=8.6;
  32.   y[3]:=14.42;    y[4]:=15.8;
  33.   y[5]:=18.92;    y[6]:=17.96;
  34.   y[7]:=12.98;    y[8]:=6.45;
  35.   y[9]:=0.27;
  36. end;        { procedure get_data }
  37.  
  38. procedure write_data;
  39. { print out the answers }
  40. var    i    : integer;
  41. begin
  42.   writeln;
  43.   writeln('  I      X       Y       YCALC');
  44.   for i:=1 to nrow do
  45.     writeln(i:3,x[i]:8:1,y[i]:9:2,y_calc[i]:9:2);
  46.   writeln; writeln(' Coefficients ');
  47.   for i:=1 to ncol do
  48.     writeln(coef[i]:8:4);
  49.   writeln;
  50.   writeln('Correlation coefficient is ',correl_coef:8:5)
  51. end;        { write_data }
  52.  
  53. procedure solve(a: ary2s;
  54.         y: arys;
  55.      var coef: arys;
  56.          nrow: integer;
  57.     var error: boolean);
  58.  
  59. var    b    : ary2s;
  60.     i,j    : integer;
  61.     det    : real;
  62.  
  63.  
  64. function deter(a: ary2s): real;
  65.     { calculate the determinant of a 3-by-3matrix }
  66. begin
  67.   deter:=a[1,1]*(a[2,2]*a[3,3]-a[3,2]*a[2,3])
  68.     -a[1,2]*(a[2,1]*a[3,3]-a[3,1]*a[2,3])
  69.     +a[1,3]*(a[2,1]*a[3,2]-a[3,1]*a[2,2])
  70. end;
  71.  
  72.  
  73.  
  74. procedure setup(var b    : ary2s;
  75.         var coef: arys;
  76.         j    : integer);
  77.  
  78. var    i    : integer;
  79.  
  80. begin    { setup }
  81.   for i:=1 to nrow do
  82.     begin
  83.       b[i,j]:=y[i];
  84.       if j>1 then b[i,j-1]:=a[i,j-1]
  85.     end;
  86.   coef[j]:=deter(b)/det
  87. end;        { setup }
  88.  
  89. begin    { procedure solve }
  90.   error:=false;
  91.   for i:=1 to nrow do
  92.     for j:=1 to nrow do
  93.       b[i,j]:=a[i,j];
  94.   det:=deter(b);
  95.   if det=0.0 then
  96.     begin
  97.       error:=true;
  98.       writeln(chr(7),'ERROR: matrix is singular')
  99.     end
  100.   else
  101.     begin
  102.       setup(b,coef,1);
  103.       setup(b,coef,2);
  104.       setup(b,coef,3)
  105.     end        { esle }
  106. end;        { procedure solve }
  107.  
  108.  
  109. procedure linfit(x,y: ary;
  110.          var y_calc: ary;
  111.          var coef:   arys;
  112.              nrow:   integer;
  113.          var ncol:   integer);
  114.  
  115. { least squares fit to a parabola }
  116. { nrow sets of x and y pair points }
  117.  
  118. var    a        : ary2s;
  119.     g        : arys;
  120.     i        : integer;
  121.     error        : boolean;
  122.  
  123.     sum_x,sum_y,sum_xy,sum_x2,
  124.     sum_y2,xi,yi,sxy,syy,
  125.     sxx,sum_x3,sum_x4,sum_2y,
  126.     denom,srs,x2    : real;
  127.  
  128. begin         { linfit }
  129.   ncol:=3;    { polynomial terms }
  130.   sum_x:=0.0;
  131.   sum_y:=0.0;
  132.   sum_xy:=0.0;
  133.   sum_x2:=0.0;
  134.   sum_y2:=0.0;
  135.   sum_x3:=0.0;
  136.   sum_x4:=0.0;
  137.   sum_2y:=0.0;
  138.   for i:=1 to nrow do
  139.     begin
  140.       xi:=x[i];
  141.       yi:=y[i];
  142.       x2:=xi*xi;
  143.       sum_x:=sum_x+xi;
  144.       sum_y:=sum_y+yi;
  145.       sum_xy:=sum_xy+xi*yi;
  146.       sum_x2:=sum_x2+x2;
  147.       sum_y2:=sum_y2+yi*yi;
  148.       sum_x3:=sum_x3+xi*x2;
  149.       sum_x4:=sum_x4+x2*x2;
  150.       sum_2y:=sum_2y+x2*yi
  151.     end;
  152.   a[1,1]:=nrow;
  153.   a[2,1]:=sum_x;    a[1,2]:=sum_x;
  154.   a[3,1]:=sum_x2;    a[1,3]:=sum_x2;
  155.   a[2,2]:=sum_x2;    a[3,2]:=sum_x3;
  156.   a[2,3]:=sum_x3;    a[3,3]:=sum_x4;
  157.   g[1]:=sum_y;
  158.   g[2]:=sum_xy;
  159.   g[3]:=sum_2y;
  160.   solve(a,g,coef,ncol,error);
  161.   srs:=0.0;
  162.   for i:=1 to nrow do
  163.     begin
  164.       y_calc[i]:=coef[1]+coef[2]*x[i]+coef[3]*sqr(x[i]);
  165.       srs:=srs+sqr(y[i]-y_calc[i])
  166.     end;
  167.   correl_coef:=sqrt(1.0-srs/(sum_y2-sqr(sum_y)/nrow))
  168. end;        { linfit }
  169.  
  170. { external procedure plot(x,y,y_calc: ary; nrow: integer);
  171. }
  172.  
  173. {$I C:PLOT.LIB }        { get ptocedure PLOT }
  174.  
  175. begin        { MAIN program }
  176.   cls;
  177.   get_data(x,y,nrow);
  178.   linfit(x,y,y_calc,coef,nrow,ncol);
  179.   write_data;
  180.   plot(x,y,y_calc,nrow)
  181. end.        { MAIN }
  182.